Formation of soliton trains in Bose-Einstein condensates 
as a nonlinear Fresnel diffraction of matter waves 

A.M. Kamchatnov 1 , A. Gammal 2 ' 3 , F.Kh. Abdullaev 4 , 
and R.A. Kraenkel 5 
1 Institute of Spectroscopy, Russian Academy of Sciences, 
Troitsk 142190, Moscow Region, Russia 
2 Instituto de Ffsica, Universidade de Sao Paulo, 
05315-970, C.P.66318 Sao Paulo, Brazil 
3 Department of Physics and Astronomy and Rice Quantum Institute, 
Rice University, Houston, Texas 77251 
4 Physical- Technical Institute, Uzbek Academy of Sciences, 
70084 Tashkent-84, G.Mavlyanov str. 2-b, Uzbekistan 
5 Instituto de Ffsica Teorica, Universidade Estadual Paulista-UNESP 
Rua Pamplona 145, 01405-900 Sao Paulo, Brazil 

February 2, 2008 

Abstract 

The problem of generation of atomic soliton trains in elongated Bose-Einstein con- 
densates is considered in framework of Whitham theory of modulations of nonlinear 
waves. Complete analytical solution is presented for the case when the initial density 
distribution has sharp enough boundaries. In this case the process of soliton train 
formation can be viewed as a nonlinear Fresnel diffraction of matter waves. The- 
oretical predictions are compared with results of numerical simulations of one- and 
three-dimensional Gross-Pitaevskii equation and with experimental data on formation 
of Bose-Einstein bright solitons in cigar-shaped traps. 

Discovery of Bose-Einstein condensate (BEC) PQ 12 E] has created new active field of 
research of quantum macroscopical behavior of matter. Among most spectacular evidences 
of such macroscopic behavior one can mention formation of interference fringes between 
two condensates [3] and creation of dark |5J E] and bright |S] solitons. The interference 
phenomenon is usually considered in framework of a linear wave theory, whereas solitons are 
treated as a nonlinear wave effect. At the same time, basically, these two phenomena have 
much in common. For example, formation of bright soliton trains in nonlinear wave systems 
is often explained as a result of modulational instability, where selection of the most unstable 
mode is a result of interplay of interference and nonlinear effects (see, e.g. HDD- Such 
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interconnection of interference and soliton phenomena is demonstrated most spectacularly 
in formation of solitons in vicinity of a sharp edge of density distribution. In this case, at 
linear stage of evolution the linear diffraction provides an initial modulation of the wave and 
further combined action of interference and nonlinear effects leads to formation of soliton 
trains. Without nonlinear effects, such kind of time evolution of a sharp wave front would 
be a temporal counterpart of usual spatial Fresnel diffraction and therefore soliton train 
formation at the sharp front of nonlinear wave can be called a nonlinear Fresnel diffraction. 

Similar formation of oscillatory structures at sharp wave front or after wave breaking in 
modulationally stable systems described by the Korteweg-de Vries equation is well known 
as a "dissipationless shock wave" (see, e.g. [ID])- Its theoretical description is given [TT] IT2*] 
in framework of Whitham theory of nonlinear wave modulations ^3] , where the oscillatory 
structure is presented as a modulated nonlinear periodic wave which parameters change little 
in one wavelength and one period. Then slow evolution of the parameters of the wave is 
governed by Whitham equations obtained by averaging of initial nonlinear wave equations 
over fast oscillations of the wave. Application of this method to modulationally unstable 
systems has been given for important particular case of soliton train formation at the sharp 
front of a long step-like initial pulse f2JE3EHEl- Here we shall consider by this method 
formation of solitons in BEC with negative scattering length (attractive interaction of atoms). 

We suppose that condensate is confined in a very elongated cigar-shaped trap whose axial 
frequency u z is much less than the radial frequency u±. In the first approximation we can 
neglect the axial trap potential and suppose that condensate is contained in a cylindrical 
trap (u z = 0) and its initial density distribution has a rectangular form. Evolution of BEC 
is governed by three-dimensional (3D) Gross-Pitaevskii (GP) equation 

h 2 1 

for the condensate wave function ip, where we use standard notation g = 47rh 2 a s /m a for 
the effective nonlinear coupling constant, a s < is the s-wave scattering length, and ip is 
normalized on the number of particles in BEC, J \ip\ 2 dr = N. For analytical treatment of 
nonlinear Fresnel diffraction it is important to determine conditions when the 3D equation 
(0) can be reduced to its one- dimensional (ID) approximation (see, e.g. [15]) 

ih% = ~2^-*** + 9id\*\ 2 *, I m 2 dz = N, (2) 

where 



g 2h a s / h 



giD — T. — 2" = a ± = \ > (3) 

Zira^ m a a ± V rn a u± 

that is the transversal degrees of freedom are frozen. It is well known (see, e.g. |10j ) 
that a homogeneous distribution with linear density uq = |\1/| 2 = const described by (J2J) 
with negative gm (a s < 0) is unstable with respect to self-modulation with increment of 
instability equal in our present notation to 

h K 

^8\a s \n -(a ± K)*, (4) 



2m a a± 
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where K is a wavenumber of small periodic modulation. The most unstable mode has the 
wavenumber 

K m ax = ^V\ a s\ n o/a± (5) 

and the corresponding increment is equal to 

Tmax = 4:\a s \n uj ± . (6) 

This means that after time ~ l/(|a s |nou;_i_) the homogeneous condensate splits into sepa- 
rate solitons (diffraction fringes) and each soliton (diffraction fringe) contains about N s ~ 
no/K max atoms. If in 3D GP equation (0) the nonlinear energy gN s K max / a\ ~ gn^/a 2 ^ in 
each solitons is much less than the kinetic energy in the transverse direction, ~ h 2 /m a a\, 
then the transverse motion is reduced to the ground state oscillations and the 3D condensate 
wave function can be factorized into tp = (f)o(x,y)^(z,t), where 0o = (v^jj -1 exp[— (x 2 + 
y 2 ) I (2dj_)] is the ground state wave function of transverse motion, and ^(z,t) obeys to the 
effective ID nonlinear Schrodinger (NLS) equation (j2J). Thus, the condition of applicability 
of ID equation (J2j) for description of solitons formation is 

n \a s \ -C 1, (7) 

which means that the instability wavelength ~ 1/K max is much greater than the transverse 
radius a± of BEC. If (JIJ) is not satisfied, then the transverse motion has to be taken into 
account which may lead to collapse of BEC inside each separate soliton. Therefore we shall 
confine ourselves to the BEC described by the ID NLS equation under supposition that the 
initial distribution satisfies the condition (|7J). 

To simplify formulae in the analytic theory, we transform to dimensionless variables 
t = 2(\a s \n ) 2 uj±t, ( = 2\a s \n z/a±, ^ = ^/2\a^\n u, so that (J2J) takes the form 

iu T + Utf + 2\u\ 2 u = 0, (8) 

and u is normalized to the effective length L of the condensate J \u\ 2 d( — L/a± measured 
in units of a±. We are interested in the process of formation of solitons (nonlinear Fresnel 
diffraction fringes) at the sharp boundary of initially rectangular distribution. Since this 
process takes place symmetrically at both sides of the rectangular distribution, we can con- 
fine ourselves to the study of only one boundary. This limitation remains correct until the 
nonlinear waves propagating inside the condensate collide in its center. If the initial distri- 
bution is long enough, this time is much greater than the time of solitons formation. Thus, 
we consider the initial distribution in the form 

*«.°>={rt 2 7>o, for c<0 « 9 > 

where 7 is the height of initial step-like distribution and a characterizes the initial homoge- 
neous phase. The problem of this kind has already been considered in some other problems 
of nonlinear physics El El QH] and we shall present here only the main results. 

Due to dispersion effects described by the second term in Eq. (JH1), the sharp front trans- 
forms into slightly modulated wave which describes usual Fresnel diffraction of atoms. In our 
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case the diffraction pattern evolves with time rather than is "projected" on the observation 
plane. The linear stage of evolution is followed by the nonlinear one in which combined 
action of dispersion and nonlinear terms yields the pattern which can be represented as a 
modulated nonlinear periodic wave or, in other words, a soliton train. We suppose that this 
soliton train contains large enough number of solitons, so that their parameters change little 
in one wavelength and one period. Then, in framework of Whitham theory, the density of 
BEC can be approximated by a modulated periodic solution of Eq. (JBJ (see [Ellin]) 

n = KC, r)| 2 = (7 + 5f - 4 7 5 sn 2 ( y/(a - (3) 2 + (7 + S) 2 6, m), (10) 

where sn(x, m) is the Jacobi elliptic function, 

= {-Vt, V = -2(a + /3), (11) 
m = 4 7 5/[(«-/5) 2 + (7 + 5) 2 ], (12) 

the parameters a and 7 are determined by the initial condition Q, and (3 and 5 are slow 
functions of ( and r. Their evolution is governed by the Whitham equation 

M +oCft4) »£+*> = 0, (13) 

where Whitham velocity v(/3, 5) is given by the expression 

V(M = -2(c + fl - w _ a)(K -j E) - - _ _ , (14) 

K = K(m) and E = E(m) being the complete elliptic integrals of the first and second kind, 
respectively. Since our initial condition @ does not contain any parameters with dimension 
of length, the parameters (3 and 5 can only depend on the self-similar variable £ = £/r. Then 
Eq. (jT3j) has the solution 

C/r = t = v(p,S) (15) 
with v(/3, 5) given by (|14j) . Separation of real and imaginary parts yields the formulae 

C/r= _ 4 ^_2(7 2 -5 2 )/( / 9-a), (16) 
(a - P) 2 + (7 - 5) 2 E(m) 



(a - /3) 2 + 7 2 - 5 2 K(m) 



(17) 



which together with Eq. ()12j) determine implicitly dependence of (3 and 5 on £ = £/r. It is 
convenient to express this dependence in parametric form 

f3(m) — a — 7 ^^(m) - (1 +mA(m)) 2 , (18) 

5(m) = 7myl(m), (19) 

where 

\ (2 - m)E(m) - 2(1 - m)K(m) , N 

v ; m 2 E(m) v ; 

Substitution of these expressions into (fTUjl . lfTTj) yields the density 77 as a function of 777. 
Since the space coordinate ( defined by Eq. (|T6*j) is also a function of m at given moment 
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r, we arrive at presentation of dependence of n on ( in parametric form. The limit m — > 
corresponds to a vanishing modulation, and this edge point moves inside the condensate 
according to the law 

(_ = (-4a + 4^27)1-. (21) 
The other edge with m — > 1 moves according to the law 

C+ = -4ar, (22) 

and corresponds to the bright solitons (or fringes of nonlinear diffraction pattern) at the 
moment r. The whole region £_ < ( < ( + describes the oscillatory pattern arising due to 
nonlinear Fresnel diffraction of the BEC with initially sharp boundary at ( = 0. 

We have performed numerical simulation of ID and 3D GP equations with the aim to 
compare approximate Whitham theory with numerical results. The ID density distribu- 
tions calculated numerically from (|SJ) and analytically are shown in Fig. 1. We see excellent 
agreement between the theoretical and numerical predictions of the height of the first soliton 
generated from initially step-like pulse, but its position given by analytical formula is slightly 
shifted with respect to numerical result. This is well-known feature of asymptotic Whitham 
approach ^T] ^] which accuracy in prediction of location of the oscillatory pattern cannot 
be much better than one wavelength. Thus, we see that the above theory reproduces the 
numerical results quite well for period of time r ~ 2. For much greater time values some 
other unstable modes different from one-phase periodic solution (|10|) can also give consider- 
able contribution into wave pattern. Nevertheless, the qualitative picture of soliton pattern 
remains the same. 

For 3D numerical simulation, the GP equation (JI} was transformed to dimensionless form 
by means of substitutions x = a±x', y = a±y', z = a±z', t = 2t'/u>±, ip' = (iV 1 / 2 /a 3 / 2 )ip, so 
that it takes the form 

fyt = — A*0 + r 2 V> - (87rJV|a a |/a_ L )|V'|V, (23) 

where primes are omitted for convenience of the notation and J \ijj\ 2 27rrdrdz = 1, r 2 = x 2 +y 2 . 
Evolution of the density distribution p(z) = J °° \if)(r, z)\ 2 2nrdr along the axial direction 
is shown in Fig. 2 for the values of the parameters corresponding to the experiment [S] 
(a s = — 3ao, u± = 2tt ■ 625 Hz, L = 300 a±) except for the number of atoms which was chosen 
to be N = 5 ■ 10 3 in order to satisfy the condition (j7J), so that \a s \n = 1.7 ■ 10~ 3 . We see that 
diffraction (soliton) pattern arises after the dimensionless time t ~ 400 which corresponds 
after appropriate scaling transformation to r ~ 2 in Fig. 1. The width of solitons in Fig. 2 
also agrees with the width predicted by ID analytical theory and numerics. The spatial 
distribution of the condensate density \ip(r, z)\ 2 is illustrated by Fig. 3. The 3D nonlinear 
interference pattern is clearly seen. For greater values of the condensate density, when ID 
theory does not apply, numerical simulation demonstrates similar evolution of the diffraction 
pattern up to the moment when collapse starts in each separate soliton. Thus, formation 
of solitons in the experiment jH] with large initial number of atoms iV ~ 10 5 goes through 
collapses with loss of atoms until the remaining atoms can form stable separate soliton-like 
condensates. The present theory emphasizes the importance of the initial stage of evolution 
with formation of the nonlinear Fresnel diffraction pattern. 
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Formation of soliton trains in BEC confined in a cigar-shaped trap has also been studied 
numerically in [T9*l I2U] . The results of ID simulation in jTU] agree qualitatively with our 
results. In numerics of j2D] strong losses were introduces to prevent fast collapse of BEC 
with large number of atoms. Nevertheless, formation of soliton trains was also observed. 

The above theory is correct for evolution time much less than period of oscillations 2-k/uj z 
in the axial trap. When the axial trap is taken into account, solitons acquire velocities in 
axial direction even if initial phase is equal to zero. The number of solitons produced 
ultimately from some finite initial BEC distribution can be found by means of quasiclassical 
method applied to an auxiliary spectral problem associated with the NLS equation JHJ) in 
framework of the inverse scattering transform method ^21 E] If the initial wave function 
is represented in the form u (() = a/^o(C) ex P(^0o(C))> then the total number of solitons is 
equal approximately to 

n s = \ J yQcH^rfC - \, (24) 

where Vq(() = d<po(()/d( is the initial velocity distribution of BEC. If there is no initial 
phase imprinted in BEC, then the total number of solitons is given by the formula 

N a = (y/2\^\/ira±) / \V\dz, (25) 

which is written in dimensional units and we have neglected a "one-half" term in ([24)1 . 

In experiment, the initial stage is usually obtained by sudden change of the sign of the 
scattering length from positive to negative one, so that initial density distribution has, for 
large enough number of atoms, the Thomas-Fermi form 

|^| 2 = (3N/4Z) (l-z 2 /Z 2 ) , (26) 

where Z is the Thomas-Fermi half-length of the condensate. Then substitution of (|26|) into 
PHjl gives 

N s = y/3N\a 3 \L/(4a ± ), (27) 

where L = 2Z is the total length of the condensate. Up to constant factor, this estimate 
coincides with one obtained in j20j by division of L by the instability wavelength 1/K max . 
Note that this estimate includes also very small solitons which cannot be observed in real 
experiments, so that it must be considered as an upper limit of the number solitons which can 
be produced from a given initial distribution. The same property of this kind of estimate for 
number of dark solitons has been observed in comparison of analytical theory with numerical 
simulations in [2Tj. The axial potential influences mainly on velocities of solitons, so the 
above estimate can be applied to the condensate in a cigar-shape trap under condition that 
inequality (J7J) is fulfilled. 

In conclusion, we have studied theoretically and numerically the process of formation of 
soliton trains near the sharp edges of the density distribution of BEC. The arising oscillatory 
regions can be considered as nonlinear Fresnel diffraction fringes of matter waves. 

This work was supported by FAPESP (Brazil) and CNPq (Brazil). A.M.K. thanks also 
RFBR (grant 01-01-00696) for partial support. A.G. would like to thank Randy Hulet and 
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Figure 1: Density distributions of BEC calculated by numerical solution of ID GP equation 
(|H|) and given by Whitham theory with initial step-like wave function Q with 7 = — 1, 
a = 0. 
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Figure 2: Density distributions of BEC p(z) along axial direction for different moments 
of time calculated by numerical solution of 3D GP equation (|23p with cylindrical initial 
distribution. 
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Figure 3: Dependence of the density distributions on radial, r, and axial, z, coordinates at 
time t = 400. 
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